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A new unsteady, cavitation model is presented wherein the phase change process (bubble 
growth/collapse) is coupled to the acoustic field in a cryogenic fluid. It predicts the number 
density and radius of bubbles in vapor clouds by tracking both the aggregate surface area 
and volume fraction of the cloud. Hence, formulations for the dynamics of individual 
bubbles (e.g. Rayleigh-Plesset equation) may be integrated within the macroscopic context of 
a dense vapor cloud i.e. a cloud that occupies a significant fraction of available volume and 
contains numerous bubbles. This formulation has been implemented within the CRUNCH 
CFD, which has a compressible “real” fluid formulation, a multi-element, unstructured grid 
framework, and has been validated extensively for liquid rocket turbopump inducers. 
Detailed unsteady simulations of a cavitating ogive in liquid nitrogen are presented where 
time-averaged mean cavity pressure and temperature depressions due to cavitation are 
compared with experimental data. The model also provides the spatial and temporal history 
of the bubble size distribution in the vapor clouds that are shed, an important physical 
parameter that is difficult to measure experimentally and is a significant advancement in the 
modeling of dense cloud cavitation. 


N UMERICAL simulation of cavitating flows and specifically the development of phase change models has 
received enormous attention from researchers. Yet all the models currently used in CFD codes have major 
drawbacks for dense cavitating flows. Cavitation models in CFD tools may broadly be classified into two 
categories: a) A bubbly framework using the Rayleigh-Plesset equation, and b) Continuum formulation solving for 
vapor phase transport. There have been numerous numerical studies using the bubbly closure model (e.g., see Ref. 
[1,2]). Here, the flowfield is seeded with bubbles and the change in bubble mass (and therefore the mixture density 
in each cell) is obtained from the Rayleigh-Plesset equations. The advantage of this formulation is that a physical 
model describes the cavitation phase change. However, the key limitation of these studies arises from their implicit 
restriction that the mixture be “dilute” i.e., the vapor void fraction ) is negligibly small. For instance, Kubota et 

al., 1 * define mixture density as p m =(1 -</> g )pi which neglects the hydrodynamics of the vapor phase. This is a 

severe restriction. In particular, for cloud cavitation problems where the vapor clouds are dense and the spatial 
extent of each distinct cloud can be large, the hydrodynamics of the vapor phase has to be resolved numerically. 

The second approach followed is a continuum formulation that is applicable to dense cavitating flow (e.g., Ref. 
[3-6]) among many other works in the literature). Here, a separate equation for the transport of the vapor phase is 
solved for and there is no restriction on the vapor volume fraction being small. The cavitation process is modeled as 
a phase change source term. The primary drawback of these prev ious w orks has been that the cavitation source term 
is specified in an ad-hoc fashion and does not account for the bubble dynamics. For instance, the barotropic model 3 
completely decouples the time scale of cavitation and enforces a user specified density-pressure relationship. The 
model by Singhal 5 does allow for finite time scale with a rate equation for phase change based on the local pressure. 
However, the key limitation here is that the neither the radius nor the number density of the bubbles in the vapor 
cloud are estimated. Since the net phase change due to cavitation results from the combination of bubble number 
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density and the individual bubble radius change, these formulations are ad-hoc and not adequate for unsteady cloud 
cavitation. 

The primary contribution of our work is the development of a fundamental model for cavitation that integrates 
the two approaches discussed above. We solve for the transport of vapor phase and are not restricted to dilute multi- 
phase flows. However, in addition we track the number density and radius of the bubbles associated in the vapor 
cloud by solving an additional equation for the total surface area associated with a bubbly cloud. This allows us to 
integrate fundamental bubble dynamic equations for changes in mean bubble radius and provides an accurate model 
for applications with extensive cavitation. The details of the formulation are given in the sections below. 

In the sections to follow, we give a brief overview of the “compressible” framework for modeling cavitating 
cryogenic fluids followed by a discussion of the new unsteady cavitation model developed. The series of rigorous 
validation tests performed to verify both the accuracy of the model to track individual bubble dynamics, as well as 
the larger macroscopic evolution of the vapor cloud are described next. The model has been tested on the following 
problems in water: 1) Unsteady bubbly pipe flow, 2) Unsteady cloud cavitation on a cavitating ogive in liquid 
nitrogen. As we discuss below, the results indicate that the model performed well and yielded fundamental insight 
into unsteady cavitation physics. 


I. Multi-Phase Equation System 

We give a very brief overview of the basic multiphase equation system here and refer the reader to Hosangadi 
and Ahuja 7 for more details. The equation system is written solved in a pressure based form as: 
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The vectors Qv, E and S are given above. 


The matrix r (= dQ / dQ v ) defines the transformation from the 


conservative to primitive variables and may further be preconditioned to obtain an efficient time-marching scheme 
(we refer the reader to Ref. [7] for more details). 

The source term for the vapor phase arises due to cavitation where m, is the net rate of vapor mass generation (or 
condensation), and the corresponding source term for the energy equation is given as m,h fg where h fg is the change 
in enthalpy resulting from the phase change and is a function of the local fluid temperature. These phase change 
source terms are discussed in a later section. 

The mixture density, enthalpy, and vapor porosity are related by the following relations locally in a given cell 
volume: 


Pm - Pgfig + P/A 

aA =pAK + pAA 

1 = <Pg + A 


(3) 

(4) 

(5) 


where p g , p L are the physical material densities, while h g and h , are the sensible enthalpy of the vapor and liquid 

phase respectively, and in general are functions of both the local temperature and pressure. In our study here, these 
properties were generated from the Standard thermodynamic database 12 available from NIST for pure fluids 8 . The 
thermodynamic properties of the fluid where specified using the saturation values from the table corresponding to 
the local temperature of the fluid. Equations (1-2) represent a stiff system with large variations in the acoustic speed 
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that are a function of the local multi-phase composition. Preconditioning techniques are used to overcome this 
stiffness and obtain an efficient numerical scheme 7 . 


II. Unsteady Cavitation Model Using Cloud Surface Area Equation 

The new, unsteady cavitation model developed incorporates formulations for bubble dynamics (e.g. Rayleigh- 
Plesset equation) as an integral part of dense cloud cavitation. By drawing on our experience in dense spray 
combustion 9 , we track the surface area associated with the cavitation cloud as it evolves spatially and temporally. 
We note that this Eulerian procedure has been formulated within the premise of a dense bubble cloud where a large 
number of bubbles are present. Bubbles comprising the cloud will be characterized by their Sauter mean radius, 
which preserves the ratio of the surface area and volume for the bubbly cloud. Thus, the source term for cavitation 
will have two independent factors controlling it: a) the net surface area given by the bubble Sauter mean radius, and 
b) the rate of change of the radius, which may be specified using the Rayleigh-Plesset equation. 

The additional equation (shown in 1-D form for discussion purposes) for the bubble surface area is given as: 
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Here s, is the source terms to the cloud surface area equation derived in a manner consistent with the Rayleigh 
Plesset equation and is defined later. For clarity, we repeat the vapor mass conservation, which was already 
included as part of the original multi-phase system (Eqn. (2)), 
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The dependent variables S g and <f> y may be written as 


S g = N * 4.0* 7r* r 2 , </> g = jV*4.0*/r*r 3 /3.0 


( 8 ) 


where N is the number density of bubbles, and r the Sauter mean radius of bubbles at each grid point. Thus, the 
Sauter mean radius, r, that conserves cloud surface area and volume is defined as, r = 2>*<j> g IS . Note that no 

restriction is placed on the vapor volume fraction being dilute. Furthermore we do not restrict the number density of 
bubbles to be a constant in each cell, and allow them to convect as governed by the hydrodynamics. 

The source terms for the net vapor mass transfer and surface area change are given by 
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where r, is the rate of change of the radius. To specify the rate of change of radius r ( , the Rayleigh-Plesset 
Equation could be specified as follows: 
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However this would require the solution of an additional transport equation for the solution of (dr/dt). In our 
study a simplified and approximate form for the radius change is used: 
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where the sign of the radius change term depends on whether the bubble is growing or decaying and is dictated by 
the sign of the P v -P term. Here, P v is the vapor pressure, while P is the instantaneous pressure at the grid node 
point at which the cavitation source term is being computed, and /?, is the liquid density. We again emphasize that 
Eqn. (11) provides the phase change for an individual bubble only while the net mass transfer is given by Eqn. (9) 
by accounting for the totai surface area of bubbles present at a given spatial grid location. Hence a more 
fundamental unsteady cavitation formulation is derived. 
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III. LES Models 

Unsteady simulations presented in this paper have been computed with a one-equation LES model unless 
otherwise noted. The one equation LES model solves for a transport equation for the sub-grid kinetic energy. This is 
then used to evaluate the eddy viscosity that is used to model the unresolved SGS stresses. The model equations are 
as follows 
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where, k is the Favre filtered subgrid kinetic energy, u, is the resolved scale velocity, p is the filtered density, 
= P(U I U J -U,Uj) is the subgrid scale tensor, v T is the subgrid scale eddy viscosity and D k is the 
dissipation of the SGS kinetic energy. The subgrid scale stresses are modeled as 
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The eddy viscosity is computed as 

v T = C v k U2 A ( 14 ) 

where, C v is the modeling constant and has been computed to be equal to 0.0667 (Arunajatesan l0 ). 

IV. Validation for Unsteady Test Cases 

Rigorous validation of the cavitation model has been performed for unsteady simulations. The following test 
cases are presented here: 1) Bubble dynamics in a fluctuating pressure environment for validating the cloud surface 
area equations; 2) Unsteady cloud cavitation on an ogive in liquid nitrogen. 

A. Bubble Dynamics in a Fluctuating Pressure Environment 

Unsteady, one-dimensional, bubbly pipe flow is analyzed; the stagnation pressure at the inlet is constant while 
the back pressure fluctuates periodically. For dilute bubble concentrations, an analytical solution is derived for the 
pressure and velocity of the fluid. The bubble response at various forcing frequencies is obtained by injecting 
bubbles at the inflow and tracking their growth/decay using a Lagrangian procedure where the pressure and velocity 
at each instant is specified using the analytical solution; this is referred to as “exact” bubble solution. The numerical 
solution is obtained using the Eulerian CFD framework. Both the mean liquid solution and the bubble dynamics are 
computed. Thus rigorous validation of bubble dynamics is achieved with spatial and temporal variations of the 
mean flow. 

The analytical solution for the velocity and pressure fluctuations is derived to be (Athavale 1 '), 
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where, Q is the non-dimensional frequency, and £ is the amplitude of the pressure fluctuations enforced at the 
exit. Note that the velocity fluctuations are uniform in space, while the pressure fluctuations show a linear decay 
from the exit to the inflow as a function of the frequency of forcing. 

The temporal history of the solution is shown in Figure 1 for three different non-dimensional frequencies of 
£2 = 1 , 10, and 100. In all these cases, the amplitude of the pressure fluctuations at the exit was kept at ten percent 
of the mean pressure. Figures la and lb show the pressure fluctuation at the inlet and exit, while Figure lc shows 
the corresponding bubble history at the exit. Excellent comparisons are obtained overall. At the lowest frequency 
( £2 = 1 ), the pressure drop between the inlet and exit is the lowest. As the frequency increases, the pressure drop 
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increases significantly and by Q = 100 , the pressure fluctuations at the inflow are almost negligible. The temporal 
history of fluctuations in bubbles radius at the exit of the pipe shows a similar behavior; with the bubble growth 
being most responsive at the lower frequencies. Indeed at the highest frequency (Q = 100), the bubble convects 
with negligible change in radius even though the amplitude of the local pressure fluctuation at the exit is unchanged. 
We caution that these results were obtained for pure vapor bubbles. For non-condensable gas bubbles, the physics is 
far more complicated since the bubbles have a natural frequency and may show resonant behavior. 


Reservoir 


P* = Constant - 


Flow 


P=P h (l+Esinwt) 





(a). Inlet Pressure History. (b). Exit Pressure History. 

Figure 1. Flow Solution History For Bubbly Unsteady Pipe Flow, £2=1, 10, 100. 


B. Cavitating Ogive Simulations 

We perform unsteady simulations of experiments by Hord 12 on a cavitating ogive to validate the formulation for 
cavitation in cryogenic fluids. Hord performed sub-scale tests using both liquid nitrogen and hydrogen in a blow- 
down tunnel. The details of the tunnel and the ogive geometry are given in Figure 2. The tunnel width is 1.158 
inches while the ogive diameter is 0.357 inches. Hence considerable blockage effects from the tunnel wall are 
present and the simulation modeled the tunnel geometry. To ensure that the tunnel blockage effects were being 
correctly modeled we performed single-phase, non-cavitating simulations and have compared it to Hord’s non- 
cavitating data in Figure 3. Good agreement is obtained giving confidence that the tunnel interaction is being 
captured. 



(a)Tunnel configuration 


(b)Ogive test article 


Figure 2. Details of Subscale Tunnel and Ogive Test Article From Hord (1973 (b)). 


Simulations are performed at the following flow conditions in liquid nitrogen: inlet temperature of 83 K, inlet 
velocity of 23.5 m/s, inlet pressure of 2.87 atm, with a freestream vapor pressure of 1.875 atm. Unsteady 
calculations for the following three different viscosity models are compared: 1) RANS k — e model, 2) one equation 
FES model, and 3) no sub-grid viscosity (or laminar) model. Fig. 4 (a-d) shows the flowfield contours for the 
unsteady calculation computed with the RANS model. As is apparent from these figures the flow remains 
essentially steady despite the calculation being performed in an unsteady mode; the vapor volume fraction and 
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temperature (Fig. 4a and 4c) show a uniform variation with the largest temperature depression occurring at the 
leading edge. The steady flowfield results from large turbulent 
viscosities in the cavity closure region (Fig. 4d) that damps the 
strength of the reentrant jet which as we shall discuss subsequently 
plays a key role in driving the cavity unsteady. The comparisons of 
pressure and temperature depression are compared with experimental 
data in Fig. 5. The experimental temperature values did not recover to 
the freestream temperature of 83°K and instead asymptoted to a value 
of 82.7°K. It was unclear if this difference was due to error in 
measurement of either the freestream or the ogive temperature, or 
whether this was caused by heat transfer losses. These values were 
shifted by 0.3°K to match the stated freestream value of 83°K. We 
observe that the results that are consistent with steady calculations; the 
leading edge temperature and pressure depression compare well with 
data but the cavity length is smaller as noted by the more abrupt Figure 3 Pressure Coefficient for 

temperature and pressure recovery. Sing|e phase Flow in Su b-Scale Tunnel. 




(c) Temperature (d) Turbulent viscosity 

Figure 4. Flowfield Contours for Cavitating Ogive Using a RANS Turbulence Model. 



with RANS Turbulence Mode! with Data, 
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(g) 


(i) 


(h) 

Figure 6. Instantaneous Vapor Volume Fraction Contours for 
Unsteady LES Calculations Shown at Intervals of 0.29ms. 

Simulations at the same conditions but computed with the one-equation LES model rather than a RANS 
turbulence model generated large scale flow unsteadiness. Figure 6 shows snapshots of vapor volume fraction at ten 
snapshots which are 0.29 ms apart (or roughly 1/1 0 th of a 340 Hz frequency). A periodic shedding is observed where 
the main cavity appears to get sheared and lifted-off the leading edge before being convected downstream. Figure 7 
shows the vapor volume fraction contours overlaid with the velocity field. It indicates that the shearing action arises 
due to fluid flowing back in a narrow region near the ogive surface. Thus the interaction of the reentrant jet with the 
boundary layer of the oncoming flow appears to play a critical role in determining the unsteady characteristics of the 
cavity. 

The instantaneous temperature and pressure contours are shown in Figures 8 and 9 respectively. One 
consequence of the unsteady shedding is that the minimum temperature is not necessarily at the leading edge but at 
the core of the traveling vapor cloud as it convects downstream. In fact the temperature at the leading edge location 
rises sharply each time the cavity gets sheared off and this will be reflected in the mean temperature profiles that we 
discuss later. The pressure contours reflect the highly dynamic flowfield where local pressure excursions are 
observed as condensation occurs in the region between vapor clouds. A detailed analysis of the frequency content 
and amplitudes of the pressure fluctuations will be provided later in this section. 

The time-averaged mean solution of the instantaneous flowfields taken over a large time sample are plotted in 
Fig. 10(a-c). Despite the large scale unsteadiness exhibited by the flow, the mean solutions look qualitatively similar 
to the steady solutions obtained with the RANS turbulence model. However some significant differences are also 
apparent. The cavity length is much longer and the variation of temperature and volume fraction much less in the 
LES mean solutions. A comparison of the temperature and pressure depression is compared to the experimental data 
as well as the RANS solution in Fig. 1 1. We note that both the temperature and pressure profile in the cavity closure 
region for the LES calculation show a more gradual recovery similar to the data and this is clearly arising from the 
unsteady effects which result in lower temperature and pressure vapor clouds convecting downstream in a periodic 
fashion. At the leading edge however, the temperature recovers more rapidly and is higher than the data over most of 
the cavity. This may also be attributed to unsteady effects since the leading edge of the cavity is not attached and 
experiences temperature increases due to condensation. 
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(g) (h) (i) 

Figure 8. Instantaneous Temperature Contours for Unsteady LES Calculation Shown at Intervals of 0.29ms. 
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(a) Vapor volume fraction (b) Pressure (c). Temperature 

Figure 10. Time Averaged Mean Flowfield Contours from Unsteady LES Calculation. 




Figure 11. Comparison of Mean Pressure and Temperature Depression from Unsteady LES Calculations 

with RANS Solution and Experimental Data. 
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To further illustrate the strong 
sensitivity of the solution to the 
interactions between the reentrant 
jet and the boundary layer, we 
compute a third solution with no 
sub-grid viscosity where the 
boundary layer is laminar. Figure 12 
compares the time-averaged mean 
temperature and pressure depression 
profiles for the LES calculation with 
unsteady laminar calculation. We 
note that the reduced viscous 
damping amplifies the differences 
we observed in our earlier Figure 12. Effect of Subgrid Viscosity Levels on Pressure and 

comparison of the RANS and LES Temperature Depression in Cavity, 

solutions. Specifically, the laminar 

solution (brown line in Fig. 13) shows an even larger mean cavity length with a slower pressure and temperature 
recovery and compares well with data in the closure region. However, the temperature and pressure rise at the 
leading edge is quicker and the comparison with experimental data gets worse here. 

The dynamic flow fluctuations arising from the unsteady cavitation physics are examined by plotting the 
pressure, temperature and volume fraction fluctuations close to the leading edge of the cavity in Fig. 13. The 



TfeMH) Tint* Ml Tim* HI 

(a) Temperature (b) Pressure (c) Volume fraction 

Figure 13. Time Traces of Flow Variables Close to the Leadings Edge of Cavity on Ogive. 



pressure fluctuations show periodic large amplitude pressure spiking where the amplitude of the peak pressure 
excursion is approximately 100-200 % of the mean pressure. 

Although the frequency between the peaks corresponds to a 
lower frequency associated with the cavity shedding, the 
pressure excursion itself occurs over a very small time scale. 

Spectral analysis of these pressure fluctuations shown in Fig. 

14 indicates that the bulk of the energy is in the shedding 
frequency mode of 400 Hz with a rapid drop-off in energy 
amplitudes at the higher overtones at 1 KHz and beyond. The 
rms value of the fluctuations is 0.49 atm for a mean pressure 
level of 1.89 atm (i.e. 25.9 % fluctuation level) which 
indicates a substantial dynamic load that can potentially cause 
damage locally. 

The temperature and volume fraction time history plotted 
in Fig. 13 indicate that when large scale instabilities occur the 
temperature fluctuations are of the order of the mean Figure 14. Frequency Spectrum of Pressure 
temperature depression itself and are driven by the Oscillations Near the Leading Edge of Cavity, 
hydrodynamic interactions between the vapor cloud and the 

pressure field. Consequently, the temperature fluctuations appear to be in-phase with the pressure fluctuations. We 
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note that the large amplitude of the temperature fluctuations results in the average temperature rising and this was 
evident on the mean temperature profiles compared in Fig. 12. The volume fraction time history confirms that 
shedding of the vapor clouds results in periods when the ogive is “wetted” by liquid nitrogen even at a location close 
to the leading edge of the cavity. 


V. Conclusion 

An unsteady cavitation model for dense cloud cavitation that can predict amplitudes and frequencies of dynamic 
pressure loads, and operate within a generalized compressible “real” fluid framework, has been demonstrated. A 
fundamental formulation for the unsteady phase change that accounts for the bubble size and number density as well 
its coupling to the acoustic propagation in the multi-phase fluid mixture has been developed. Simulations of a 
cavitating ogive in liquid nitrogen indicate that unsteady shedding of cavities appears to be controlled by the 
interaction of the reentrant jet with the vapor cavity. Highly dynamic flowfields with large dynamic pressure 
fluctuations were observed with the appropriate treatment of the turbulence terms. The validity of the unsteady 
solutions was verified by time-averaging the solution and comparing mean profiles to experimental data; the 
comparison in the cavity closure regions improved significantly within an unsteady framework i.e. a combination of 
the new unsteady cavitation model and the LES turbulence model. However, we note that resolving unsteady 
interactions and separation for boundary layers in very high Reynolds number flows is extremely difficult; in 
particular the leading edge of the cavity may be more steady/stable than indicated by the LES calculations here and 
may require the use of a hybrid LES-RANS model. 
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